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ABSTRACT 



This thesis analyzes the problem of motion stability of submarines in free positive 
buoyancy ascent under casualty conditions such as control surface jam and loss of 
propulsion system response. We employ fully nonlinear, coupled six degree of freedom 
equations of motion and we allow response to occur in combined vertical and horizontal 
planes. Continuation and homotopy theory techniques are utilized to trace all possible 
steady state solutions in six degrees of freedom, while local perturbation reveals their 
stability properties. Vehicle geometric properties and control surface deflections are used 
as primary bifurcation parameters. Regions in parameter spaces are identified where 
extreme sensitivity of solutions to geometric properties and hydrodynamic modeling is 



present. 
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I. INTRODUCTION 



The dynamic response of a submarine under casualty conditions constitutes a 
crucial, and frequently limiting, factor in establishing the vehicle’s submerged operating 
envelope. As basic casualty conditions in the context of this work, we mean the loss of 
control surface and/or propulsion system response, and a flooding casualty where the 
boat must either be brought to the surface or stabilized to a new operating depth. Of 
particular significance is the study of the ability of the boat to recover from a control 
surface jam. There exist several factors which determine the severity of such a situation 
as well as the recovery procedures; the initial conditions (speed, depth, pitch angle, etc.) 
during the jam, the actual control surface angles, reversing time and backing power of 
the propulsion system, the ability to blow ballast, and the time between recognition of 
casualty and initiation of proper recovery procedures. To this end, it is crucial that we 
have a clear understanding of the dynamics of the boat during an emergency situation. 

The traditional methods for establishing dynamic stability of motion concentrate 
mainly on eigenvalue analysis during small perturbations around nominal straight line 
paths Clayton and Bishop [Ref. 1]. Two indices are utilized, a stability index Gy for the 
vertical plane and G|, for the horizontal plane Roddy [Ref. 2]. In terms of the slow 
motion derivatives, these indices are given by: 
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Z^(mg-XQm) 

* YjN^-x^m) 

Positive values for these indices, whose usage dates back to the 1950’s, indicates motion 
stability in the corresponding plane. The underlying assumption in these criteria is that 
during normal straight line motions, the coupling between horizontal and vertical plane 
motions is relatively weak, and can be neglected. Although vortex shedding and flow 
separation introduces a certain degree of coupling, the above assumption has been proven 
quite useful in design and analysis. However, for a high-speed fast-maneuvering 
submarine operating at the extremes of her submerged operating envelope, or during 
emergency situations, the above assumption of uncoupled motions breaks down. High 
amplitude motions may take place in all six degrees of freedom, and the nonlinear 
interactions between the various modes of motion become more pronounced. Therefore, 
we have to carefully consider the motion characteristics allowing for coupling between 
horizontal and vertical planes. 

The implications of nonlinear effects and coupling are numerous. In the case of roll 
motion, which is one of the most critical responses, there is growing evidence of 
complicated dynamics and chaotic response under certain excitations Falzarano, Shaw and 
Troesch [Ref. 3], Taz U1 Mulk and Falzarano [Ref. 4]. 
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The motion of a submarine in the ocean environment involves some of the most 



complex fluid-structure interaction problems. For example, the steady and unsteady 
characteristics of the stem, and especially with regard to the effect of the wake and 
rolled-up body vortices on propulsor unsteadiness, and consequences for noise and 
hull vibrations are important problems of great concern. At present, too little is 
known in depth about the flow at the stem of ships and submarines and the 
footprints of their ensuing wakes in homogenous and stratified media Sarpkaya 
[Ref. 5]. 

Typical vortex stmcture about a submerged body can be seen in Figure 1 , which is taken 

from Lugt [Ref. 6]. 




Figure 1. Typical Vortex Stmcture About a Submerged Body (From Lugt, 1981). 

When an axisymmetric body moves at a sufficiently large angle of attack, the 
boundary layer vorticity may lift olf the surface and sheets of vorticity are 
convected downstream while rolling up into streamwise vortices on the leeward 
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side of the body. This separation and roll-up phenomenon, known as the crossplane 
separation, occurs on almost all maneuvering bodies at sufficiently high angles of 
attack. The roll-up of the vortices and their subsequent evolution are extremely 
important for the proper design of a submerged body and more important for the 
determination of the interaction of the vortices with the stem and the propulsion 
system. It is this interaction that determines the maneuverability and control of a 
submerged body Saipkaya [Ref. 5]. 

In the case of submarine motions, there is evidence of bifurcation phenomena and 
extreme sensitivity of response to initial conditions and control actions during emergency 



TYPICAL SIMULATION 




Figure 2. Typical Simulation Results. 
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ascent scenarios such as recovery from a dive plane jam. As motivation for the analysis 
that follows we present typical simulation results in Figure 2. The time simulation is in 
terms of vehicle roll angle versus time for 2 % excess buoyancy with the buoyancy 
force located 1 % of the vehicle length forward of the center of gravity, -6 degrees dive 
plane angle, and 0. 1 feet metacentric height. Four different recovery actions are shown, 
all parametrized by the applied rudder angle in degrees. It can be seen that although 
zero rudder angle appears to bring the vehicle roll angle back to zero in the shortest 
time, the response does not persist. Small nonzero values for the rudder angle develop 
excessive roll angles in a divergent way, while larger rudder angles reduce the amount 
of roll. Clearly the vehicle effective rudder angle is largely affected by the amount of 
vorticity and currents in the flow field, and its actual value can not be known exactly. 
Since situations like the one presented in the graph should be avoided, we need to 
develop a mechanism for assessing those regions in the parameter space where the 
response can not be simulated with confidence. 

Stability in emergency ascent has been studied in the vertical plane. In Booth [Ref. 
7] and [Ref. 8], the vehicle response was distinguished into either a nearly vertical ascent 
or predominantly forward motion. In Papoulias and McKinley [Ref. 9], it was found that 
the above distinction is not always meaningful as a result of the many parameters that 
affect the problem. This latter study maintained some vertical plane restrictions, although 
it indicated the potential existence of pitchfork bifurcations which led to coupled out-of- 
plane solutions. In this work, we relax the requirement for vertical plane motions, and 
we analyze the stability properties of all possible steady states in six degrees of freedom. 
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We employ a combination of singularity theory Golubitsky and Schaeffer [Ref. 10], 
bifurcation theory Guckenheimer and Holmes [Ref. 11], and numerical continuation 
methods Seydel [Ref. 12] in order to capture all steady state solutions that are physically 
admitted by the coupled nonlinear equations of motion. The primary bifurcation 
parameters used are: amount and location of excess buoyancy, diveplane and rudder 
deflection, and the metacentric height. Solution branching is shown to occur in various 
forms, including single and multiply connected pitchfork bifurcations, separation of 
solutions, hysteresis, and teardrop branches. Dynamic loss of stability in the form of 
Poincare- Andronov-Hopf bifurcations to periodic solutions is also identified. It is shown 
that for certain ranges of parameters, these periodic solutions persist in the vicinity of 
inverted pendulum steady states. Finally, we summarize our results in the form of 
bifurcation graphs which identify parameter regions with qualitatively different 
asymptotic response characteristics. For demonstration purposes, all computations in this 
work are performed for the Swimmer Delivery Vehicle, a 17.4 feet vehicle, for which 
a complete set of hydrodynamic and geometric properties is available Smith, Crane and 
Summey [Ref. 13]. Unless otherwise specified, all results in this work are presented in 
dimensional form, linear dimensions in feet, velocities in ft/sec, angular deflections in 
degrees, and time in seconds. 
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n. PROBLEM FORMULATION 



A. PROBLEM STATEMENT 

Controlling emergency ascent situations on submersible vehicles such as dive plane 
jam recovery is of concern to the submariners. In order to control such situations, one 
must first be able to predict the dynamic response of positively buoyant submersibles. 

Dynamic response equations of motion describe the maneuvering characteristics of 
submersible vehicles for six degrees of freedom. These equations assume constant 
coefficients for hydrodynamic forces and moments approximated by zero frequency added 
mass and damping terms plus the quadratic terms for drag forces. The constant 
coefficients vary for each vehicle and depend on such things as vehicle body shape, 
location and magnitude of vehicle buoyancy, position of bow and stem planes, position 
of rudder, vehicle speed, vehicle mass characteristics, vehicle hydrodynamic coefficients, 
propeller ipm and control surface inputs. This thesis uses the equations of motion and 
hydrodynamic coefficients for a submerged Mark IX Swimmer Delivery vehicle (SDV) 
developed by Smith, Crane and Summey [Ref. 13] to forecast the dynamic behavior of 
submersibles in a positively buoyant condition. 

For submarines where high amplitude motions may take place in all six degrees of 
freedom, nonlinear interactions between the various modes of motion may become more 
pronounced. In particular, there is growing evidence of bifurcation phenomena during 
high speed maneuvering and emergency ascent scenarios such as recovery from a dive 
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plane jam. In these cases, it is no longer true that decoupled linear analysis techniques 
are sufficient and one is forced to consider the tme character of six degrees of freedom 
motions. So, in this thesis, six degrees of freedom equations of motions have been used 
to investigate branching behaviors of the equations of motions. This is done under the 
condition of changing certain system parameters while keeping others constant. 

B. EQUATIONS OF MOTION 

The equations of motion given below are referenced to a right-hand orthogonal axis 
system fixed in the body as shown in Figure 3. The origin of this body axis system is 
located in the vehicle’s X-Z plane of symmetry on the mid-body centerline behind the 
vehicle’s nose. Positive directions for control surface deflections, forces and angular 
moments, and linear and angular velocity components are shown in Figure 3. In the 
following equations, differentiation with respect to time is denoted by a dot over a 
quantity. The six degrees of freedom equations of motion for a submarine in surge, 
sway,heave, roll, pitch, and yaw, respectively, are: 

mx [u-vr+wq-X(;{g^+r^) +y^{pg-t) +z^{pr+^) ] =Xff+X„+Xc (D 
mx [v+ur-wp+Xf,(pg+r) ^z^igr-p) ] =Y^+y',,+y<, (2) 

mx [w-ug+vp+X(,ipr-g) +yc(gr+p) -z^{p^+g^) ] =Zff+Z„+Zc O) 

J^p+ (I^Iy) gr+J^(pr-g) -ly^(g^-z^) -I^^{pg+r) 

+mx [y^iw-ug+vp) -z^iv+ur-wp) ] =K„+K,f+K(, 
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roi* viKw 
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Figure 3. Schematic Indicating Positive Directions of Axes, Angles, Velocities, Forces, and Moments 



Iyq+{I^-I^)pi-I^{qz+p) +Iy^(pq-r) +I^^{p^-r^) 

-mx [x^iw-uq+vp) -z^iu-vz+wq) ] =M^+Af^+Af(, 



( 5 ) 



-Iry(P^-Q -IJpr +4) +IJqr-p) 

+mx[Xf/y+ur-wp) -yff/u -vr+wq)] =N„+Njy+N^ 

In these equations, the left hand sides represent inertial forces and moments 
(Newton’s Law) and the right hand sides model the external forces. Subscript H reflects 
hydrodynamic contributions, W buoyancy and weight effects, C forces arising from 
control surface (rudders, dive planes, and bow planes) actions, and the rest of the 
symbols are based on standard notation and will be explained at the end of this section. 
The hydrodynamic radiation and viscous forces are expressed as: 






YH=YpP*Y/*yp^q^ygr^r*Y.v*Y^up*Y^ur+Y^vq*Y^wp*Y^wr-*-Y^uv*Y^^ 

Z„=Z^q+Zpy+Zp^r+Zy^ZjiV+Z^uq-^Z^vp+Z^vr+Z^uw^Z^v^ 

-|p A(x)(v+xr) 2+C^^Mx)(w-x9)2]i^^^ 



K„=K^*K/+Kp^q +K^/ir^K.v+Kjfip ^K^ur+K^vq 
+K^wp*K^wr+K^uv*K^vw 



(10) 
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( 11 ) 



*^piP ^ ^ *M^uq +M^vp +M^vr +M^uw +M^ 

+ ^ P f^'^[Cj,h(x)(v+xrf +Cj,JKx)(w-xq)^]^^^xdx 



-|p ly[Cj,h(xKv^xrf^Cj,^Kx)(w-^^^^ 



These are given in customary form of series expansion in terms of the 
hydrodynamic coefficients and cross flow integral terms which are integrated over the 
entire length of the body and represent quadratic forces. The cross flow velocity U,f is: 

U^=[(v+xrf-^(w-xqfy'^ (13) 

Hydrodynamic restoring forces and moments are due to the vehicle weight W and 
buoyancy B, and are given by: 



X^=-iW-B)sinB (14) 

y^=()f-B)cosesin<|> (15) 

Z^=(ff-B)cos0cos<l> (16) 

^ij,=(yg)r-y^)cos0cos4i -(z^lT-ZgB)cos0siii4> (17) 
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M^= -(x^lF-a:^)cos6cos4> -(Z(jlF-z^)sm0 



(18) 



N^={XfjW-x^)co&Q^'^ +(y^PK-y^)sin0 (19) 



Forces and moments, due to control surface deflections, are reflected as added drag 
in surge, while in sway, heave, pitch, and yaw they are directly proportional to the 
control surface deflections. 






( 20 ) 



( 21 ) 

Zc=u\Z,b^^Z,6^ (22) 

Kc=0 (23) 

Mc=u^(A/,6,-hM,6^ (24) 

Nc=N^u^6^ (25) 



Usually, control surface deflections are kept intentionally small, and the linearity 
assumption in Equations (21), (22), (24) , and (25) remains valid. Unlike the surface ship 
case, the roll moment Kc is zero for a submersible since the rudder is centered with the 
vehicle center plane. The hydrodynamic coefficients in equations (7) through (12) are 
functions of the frequency of motion, or what amounts to the same thing functions of 
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the maneuver at hand. In this work slowly varying reference motions are studied and 
therefore, it is assumed that they remain constant and equal to their zero frequency limit. 
It should be emphasized though, that the constant coefficient assumption would break 
down in studies related to the fast motions under the action of first order wave forces in 
the case of a nearly surfaced submarine. Another important assumption in this study is 
that they are assumed to be constant throughout the range of vehicle angle of attack. 
Ordinary maneuvering models are usually validated for angles of attack + 15 degrees. 
For higher angles of attack, the cross flow drag terms CpY and Cpz dominate the 
response, and they are functions of the side and angle of attack. Considering them to be 
constant does not alter significantly the behavioral characteristics and qualitative 
bifurcation results that are derived. Similarly , Cpy and Cp^ are functions of speed due 
to the Reynolds number effect on cross flow drag. This is more pronounced in small size 
unmanned untethered vehicles, whereas for submarines the cross flow drag terms remain 
relatively constant over the entire speed range. 

Since the equations of motions are referred to an axis system that is fixed in the 
vehicle, and thus translates and rotates with it, the orientation and the position of the 
moving body axis system relative to fixed inertial reference system must be specified. 
The orientation of the body axis system with respect to the inertial reference system is 
defined by the standard Euler Angles ’i'(yaw), ©(pitch), and #(roll). The rotation 
sequence from the inertial reference system to the body axis system is 0, $ as shown 
in Figure 4. So to complete the model, we need the expressions for the Euler Angle 
rates of change; 
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(1) Vehicle-Centered 
Gravity-Directed System 
parallel to inertial 
axis system. 



(2) Horizontal Flight 
Reference System derived 
from X(|YoZo by rotation a- 
bout Z through yaw angle ii. 




from X*'Y*'Z“ by rotation 
about Y" through pitch 
angle 9. 



erence System derived from 
X'Y'Z' by rotation about X' 
through roll angle a. 



Figure 4. Unit Sphere Development of Euler Angles 
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4)= p+ q sin<|) tan6 + r cos<() tan0 



(26) 



Q=q cos<t)-r sin<|) 



(27) 



,jr= q 

COS0 COS0 



(28) 



Finally, it is assumed that propulsion is inoperative and the propeller is rotating 
freely. Therefore, propulsive forces are not included in Equations (1) through (6). The 
driving mechanism for the vehicle is its excess buoyancy, B-W > 0. The problem is then 
to assess the asymptotic dynamic characteristics of the system during this condition of 
free positive buoyancy ascent. 

The model presented above can be written in its state space form by selecting as 
state variables; 



where the first six describe the system motion and the last two its geometry. Notice that 
the yaw angle does not affect the equations of motion and is, therefore, not included 
in (29). Angle can be computed from (28), once the time histories of the state 
variables have been obtained. 

In a compact vector form the state equations can be written as. 



x^=u , , x^=w , x^=p 

Xs=<l » x^=r , x,=<|) . Xj=0 



(29) 



X = /(X) 



(30) 



where bold face indicates a vector. 
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Major variables and parameters as defined by Smith, Crane, and Summey 



[Ref. 13] are given below: 

1. Dynamic Variables 



u,v,w 


- Linear velocity components of vehicle with 
respect to orgin of body axes system relative to fluid. 


p,q,r 


-Angular velocity components of vehicle with 
respect to body axes system relative to inertial 
reference system. 


X,Y,Z 


- Hydrodynamic force components along body 
axes. 


K,M,N 


-Hydrodynamic moment components along body 
axes. 



2. Mass Distribution Parameters 



m 


- Mass of the flooded vehicle, including the mass of the 




entrained fluid. 


W 


- Weight of the flooded vehicle , including the weight of 
the entrained fluid. (=gm ;where g is the 
acceleration of gravity). 


V 


- Displacement volume of the vehicle. 


B 


- Buoyancy force acting on the vehicle (pgV). This is 
independent of the inertial mass distribution of the 
submersible vehicle, including whether or not it is 
flooded . 


> Yq , Zq 


- Coordinates of the CG (center of gravity) in the body 
axis system (Figure 3). These will depend on the mass 
distribution of the vehicle, including the mass of the 
entrained fluid. 


Xb, Yb, Zb 


- Coordinates of the (center of buoyancy) in the body axis 
system (Figure 3). These are independent of the mass 
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distribution system, but may vary with the addition or 
removal of external appendages. 

I„ ly, - Moments of inertia about the body system axes, 

including the entrained fluid. 

lyz • Products of incitia about the body system axes, 
including the entrained fluid. 

3. Remaining Parameters 

- Mass density of fluid medium. 

- Width and height of vehicle in its xy and xz 
planes, respectively, at location x measured in the 
body axes system (Figure 3). These quantities are 
required in the integral defining the cross flow 
forces and moments in the equations of motion. 

- Coordinates of vehicle nose and tail as 
measured in body axis system. (Figure 3). 

- Stemplane, bowplane and rudder deflection 
angles in radians (Figure 3). 

C. STEADY STATE CONDITIONS 

Steady state conditions are achieved when the submersible reaches constant linear 
and angular velocities which must assume finite values. Therefore, the body fixed linear 
accelerations and body fixed angular accelerations will be zero. Likewise, the vehicle 
will have reached constant angles of roll <l> and pitch 6 , making the derivatives of them 
equal to zero. When these values are put into Equations (1) through (6), and substituting 

<j) = 0 = 0 

in Equations (26) through (28) the steady state values of the angular velocities can be 
found as: 



p 

b(x), h(x) 



^o*e» ^^»il 

5s AA 
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p = -iji sin0 


(31) 


q = i|r sin<t> cos6 


(32) 


r = \|f cos(|» cos e 


(33) 



Combining Equations (31), (32), (33) with the equations (1) through (6) yields a 
system of nine coupled nonlinear algebraic equations in the form 

fix) =0 (34) 

where the overbar denotes an equilibrium solution. 

Equation system (34) can be solved for the steady state values of u, v, w, p, g, r, 
4>, 6,'^. This is a highly nonlinear system of equations, and it may exhibit solution 
branching and/or multiple solutions McKinley [Ref. 14]. Here ^ is not zero at steady 
state because taking ^=0 at steady state will restrict the analysis to the vertical plane. 
In this analysis, motion in all six degrees of freedom is studied. 

Equation system (34) can be written in the following form 

/(;,i)>0 (35) 

Here, X denotes any one of the system parameters. By system parameter, we mean the 
effects of control surfaces, propeller rpm the values of buoyancy and weight, etc. 
Therefore, our system of Equations (35) represent a multiparameter problem. This 
multiparameter problem can be solved by keeping all parameters fixed except one 
(denoted by X). This system of equations for various values of X can be solved 
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continuation method. The continuation method will be explained in the foUowing 
section. 

D. PRINCIPLES OF CONTINUATION 

The system of nonlinear ’algebraic’ equations 

/ ( y , A ) = 0 (36) 

serves as a basis for the discussion. Here y denotes an n-dimensional vector. No 
generality is lost by restricting attention to Equation (36). Both steady state solutions of 
ODE’S and PDE’s are solved by approximating them by such system of nonlinear 
equations. 

Assume that at least one solution of Equation (36) has been calculated. Let us 
denote this first solution by (y* , X|).The continuation problem is to calculate further 
solutions, 

. ^2) > (y^ » ^3) » • • • • 

until one reaches a target point, say X = X^. (The superscripts are not to be confused 
with exponents.) 

The j-th continuation step starts from (an approximation of) a solution (y*, Xj) of 
equation (36) and attempts to calculate the solution (y"^', Xj+,) for the next X, namely 

^j+i* 

With predictor-corrector methods, the step j-*j + l is split into two steps: 
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Predictor step Corrector step 

see Figure 3 Seydel [Ref. 12]. 




Fig:ure 5. Schematic Showing Predictor Corrector Method (From Seydel, 1988) 

In this, thesis a predictor corrector method is used. In general, the predictor is not 
a solution of equation (36). The predictor merely provides an initial guess for corrector 
iterations that home in on a solution of equation (36). The major portion of the work is 
either on the predictor step (resulting in an approximation close to the branch) or on the 
corrector step (if the predictor has produced a guess far from the branch). The distance 
between the two consecutive solutions (y^ Xj) and (y*^*, Xj+,) is called the step length or 
step size Seydel [Ref. 12]. 

It is obvious that first of all, we need a solution procedure for the particular class 
of problems under investigation. In this work, a Newton-like method has been used to 
find the solutions of the problem. Therefore, calculation of the first solution is often the 
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most difficult part. Initial guesses may be based on solutions to simplified versions of 
the underlying equations. Sometimes choosing special parameters (often equal to zero) 
leads to equations that are easily solved. In complicated cases, a hierarchy of 
simplifications must be gradually relaxed step by step, the solution of each equation 
serving as an initial guess to the following more complicated equation. Solving a 
sequence of equations with diminishing degree of simplification until finally the full 
equation solved is called homotopy Seydel [Ref. 12]. 

Equation system (35) can be simplified by assuming the angle of yaw (^), its 
derivative, and p, q, r to be zero. By doing this, the motion of the vehicle will be 
restricted to the vertical plane. In this case, the resulting equations can be solved 
analytically. These analytical solutions were done in McKinley [Ref. 14]. In this work 
these solutions and certain numerical simulation results have been used as initial 
conditions to the continuation runs. 

The continuation approach does not guarantee that all possible solutions can be 
located in a specific example. In particular, changes to detect isolated branches are 
small. Therefore, in this work there might be other solutions which had not been located 
by the continuation algorithm which is used. In this thesis, a continuation method is 
implemented BIFPACK [Ref. 15]. A careful combination of analytical techniques, 
continuation methods, and numerical simulations are used to ensure that all steady state 
solutions are captured. 
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m. CONTINUATION IN 6s 



A. GENERAL 

In this chapter, the results of continuation in terms of the stemplane angle will be 
discussed. As was mentioned before, it is very important to have an initial condition 
before starting a continuation study. In his master’s thesis, McKinley solved the system 
of equations in the vertical plane [Ref. 14]. In that work, he found that when the 
longitudinal center of buoyancy (Xqb) was equal to -1 percent of the vehicle length, for 
a certain value of stemplane, deflection the real part of one of the system eigenvalues 
became zero while the others stayed negative. By using this fact, it was assumed that 
there is a pitchfork bifurcation in our system for the above specific values of system 
parameters. In this work, the results of his thesis were used as initial conditions in the 
continuation mns to find the vertical plane solutions. Since these results were the 
analytical solutions of the simplified system equations by restricting the motion in vertical 
plane, we were able in this work to validate the continuation mns. In continuation mns, 
the system of equations Equation (35) was solved for changing values of the primary 
parameter X (in this case X=6s). After a few continuation mns, a pitchfork bifurcation 
was found. This pitchfork bifurcation is explained in the following section. 
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B. PITCHFORK BIFURCATION 



Suppose we have the nonlinear system of state equations, 

X = f (x) ( 38 ) 

we know that the equilibrium points, x of the system are defmed by: 

f{x)=0 ( 39 ) 

This is a nonlinear system of algebraic equations, and it may have multiple 
solutions in x, which means that the nonlinear system may have more than one position 
of static equilibrium. If we pick one equilibrium x , we can establish its stability 
properties by linearization. The linearized system becomes 

x=Ax , 

where A is the Jacobian Matrix of f(x) evaluated at x. 




and the state x has been redefmed to designate small deviations from the equilibrium x, 

x-x-x. 

As long as all eigenvalues of A have negative real parts, we know that the linear system 
will be stable. This means that the equilibrium x will be stable for the nonlinear system 
as well. 

The question we ask ourselves next is, what happens if one real eigenvalue of the 
linearized matrix A is zero? The interesting case here is when the rest of the eigenvalues 
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have all negative real parts, otherwise x is unstable and the problem is solved. If the 
case of a zero eigenvalue appears to be too specialized to be of any practical use, 
consider this: Assume that f(x) depends on one physical parameter, and that physical 
parameter is allowed to vary over some range. Then it is clear that A will depend on 
that parameter and as the parameter varies, it is possible that one real eigenvalue of A 
will become zero for a specific value of the parameter. Our problem is then to establish 
the dynamics of the nonlinear system as one real eigenvalue of A crosses zero; i.e. , goes 
from negative to positive. As the solutions evolve with time, things are interesting only 
along the direction of the eigenvector that corresponds to the critical eigenvalue (the one 
that crosses zero). Along the rest of the directions in the state space, everything should 
converge back to equilibrium; remember that we assumed that all remaining eigenvalues 
of A have negative real parts. 

We can see then that it is possible to approximate our original system by a one 
dimensional system, which is much easier to analyze. The dynamics of the two systems 
will be qualitatively similar. The formalization of the above reduction procedure 
constitutes what is known as center manifold reduction, or normal form computation in 
nonlinear analysis. So, let us see what happens for the case of a zero eigenvalue by 
using a (typical) first order system, 

x=Xx-x^ , 

for X > 0 there are two nontrivial equilibria, x= ± VX. The transition of stability is 
illustrated by Figure 6. We can summarize the stability as follows; 
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• For X < 0 only the trivial equilibrium exists and is stable. 

• For X > 0 the trivial equilibrium becomes unstable and a pair of symmetric stable 
equilibria are generated. 




Figure 6. Supercritical Pitchfork Bifurcation 



This phenomenon, the loss of stability of an equilibrium and the generation of additional 
equilibrium states, is called a pitchfork bifurcation and is very common in nature; Euler 
buckling of a beam is a very typical example. In particular, the above case is referred 
as the supercritical pitchfork, this is a rather benign loss of stability since upon loss of 
stability of the trivial equilibrium the additional nearby equilibrium states are stable. 
Occasionally, the above case is referred to as a soft loss of stability, since for small 
values of X beyond its critical value, the final steady state of the system does not differ 
much from the nominal (trivial) steady state. 
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As a second example, consider a "similar" system as before, the linear part remains 
the same, and the nonlinear part changes sign, 

x=Xx+x'^ 

equilibria and stability behavior for this equation are shown in Figure 7. There is again 
a loss of stability at the bifurcation point (y, X) = (0, 0), but in contrast to the first 
example of Figure 6 there is no exchange of stability. Instead, the stability is lost locally 
at the bifurcation point. We can summarize the stability as foUows : 

• For X > 0 only the trivial equilibrium exists and is unstable. 

• For X < 0 the trivial equilibrium becomes stable, and a pair of symmetric 
unstable equilibria are generated. 




Figure 7 . Subcritical Pitchfork Bifurcation 

This case, which is also shown in the figure, is called a subcritical pitchfork. A 
comparison with the first case reveals that this is a much more serious loss of stability 
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case. Upon loss of stability of the trivial equilibrium position, there is no other stable 
equilibrium. 

Now let us pertuib the second example, considering the equation 

Xx+x:3+p=0 (40) 

An analysis reveals that in the perturbed case the bifurcation is destroyed. 
Differentiating equation (40) with respect to X gives: 



V+ 1. 



dx 






Since x=0 is not a solution of equation (40), there is no horizontal tangent 
dx/dX=0 to the solution curve x(X). Checking for a vertical tangent means considering 
the dependence of X on x and differentiating with respect to x. This yields 

4^=0 for X=-3x2 , 



which, substituted equation (40), in turn yields the loci 

(x,X) = [(-|)^ ,-3(1)^] 

of the vertical tangents. There are three solutions of equation (40) for 

X ^ -3(f)^ 

and the branching diagram looks like Figure 8 (solid lines). The point with a vertical 
tangent is a turning point. For < 0, the curves in Figure 8 are reflected in the X-axis. 
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As )S -* 0, the curves approach the pitchfork (dashed curve). For an arbitrary 
perturbation the occurrence of a turning point is typical. Bifurcation is the rare 
exception because it occurs only for j8=0 . 




In our case, we found a pitchfoik bifurcation. We found that the trivial solutions 
of our system come from in-plane (vertical plane) solutions. After the pitchfork 
bifurcation, we had a branching in the steady state solutions. The physical meaning of 
these branches is that they correspond to out of plane motion of the vehicle at the steady 
state. To get these branches, we used the numerical simulation results as an initial guess 
for the continuation run. Figure 9 shows a typical pitchfork bifurcation that we found is 
a supercritical pitchfork bifurcation. In this and all subsequent runs in this chapter, 5s is 
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our main bifurcation parameter (X) and dr is the perturbation parameter (/8). In Figure 
9 and aU subsequent figures, solid lines correspond to the symmetric system (6r=0 
degrees) dashed lines to 6r=0.5 degrees, dotted lines to dr = 2.5 degrees, and dash-dot 
lines to dr= 5 degrees. As we can see in this figure, it illustrates all the typical 
characteristics of a pitchfork bifurcation. 




-20 -15 -10 -5 0 5 10 15 20 

STERN PLANE ANGLE 



Figure 9. Pitchfork Bifurcation From Continuation Runs 

C. CONDITIONS 

1. Defuimg Additional Terms 
a. Excess Buoyancy SB 

Excess buoyancy is defined as dB = B - W where B is the submersible’s 
total buoyancy and W is the submersible’s total weight 



29 



h. Longitudinal Center of Buoyancy , Xgb 

The longitudinal center of buoyancy is defined as Xgb = Xq - Xg where 
Xq is the longitudinal center of gravity with respect to the body fixed axis and Xg is the 
longitudinal center of buoyancy with respect to the body fixed axis. 

c. Vertical Center of Buoyancy, Zcb 

The vertical center of buoyancy is defined as Zgb = Zg - Zb where Zq is 
the vertical center of gravity, with respect to the body fixed axis, and Zg is the vertical 
center of buoyancy with respect to the body fixed axis. 

2. Assumed Conditions 

a. Lateral Centers of Gravity, y^, and Buoyancy, y^ 

The lateral center of gravity and center of buoyancy are assumed to be 
on the same plane(yo = ys = 0)- 

b. Propeller Speed , n (revolution per minute) 

The propeller speed is assumed to be zero (n=0). 

c. Propeller Coefficients , and 

From Smith, Crane, and Summey [Ref. 13], the propeller coefficients are 

zero(Kp,op = = 0). 

d. Vertical Center of Buoyancy, Zcb 

The vertical center of buoyancy is assumed to be positive. 
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D. STEADY STATE RESULTS 



Figures 10 through 18 show steady state solutions for the eight state variables, 
namely; u, v, w, p, q, r, 4>,0. In Figures 10 through 18, solid lines correspond to 5r = 
0 degrees, dashed lines correspond to 5r =0.5 degree, dotted lines correspond to 5r = 
2.5 degrees, and dash-dot lines correspond to 6r = 5 degrees. Figures 16 and 17 show 
steady state solutions for roll angle in two different axis scales. In all these cases, 
stemplane action is continuation parameter (i.e. ,X) and 6r is the perturbation parameter 
(i.e. ,/3). Values of the other parameters are kept fixed:excess buoyancy 5B = 2 % of 
the vehicle weight (W); deflection of bow planes, Sb = 0; location of horizontal and 
vertical centers of buoyancy, Xg = Zg = 0; location of vertical center of gravity ,Zob = 
0.1 feet;and location of longitudinal center of buoyancy; Xgb= -1 % o the vehicle length. 



SURGE VELOCITY U vs DS 




STERN PUXNE ANGLE 

Figure 10. Steady State Surge Velocity 
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SWAY velocity V vs OS 




Figure 11. Steady State Sway Velocity 



HEAVE VELOCITY W vs DS 




STERN PLANE angle 

Figure 12. Steady State Heave Velocity 
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-20 -15 -10 -5 0 5 10 15 20 

STERN PLANE ANGLE 

Figure 13. Steady State Roll Velocity 




-20 -15 -10 -5 0 5 10 15 20 

STERN PLANE ANGLE 

Figure 14. Steady State Pitch Angular Velocity 
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YAW ANGULAR VELOCITY r vs OS 




-20 -15 -10 -5 0 5 10 15 20 



stern plane angle 

Figure 15. Steady State Yaw Angular Velocity 



ROLL ANGLE PHI vs DS 




STERN PLANE ANGLE 

Figure 16. Steady State Roll Angle 
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THETA 




-20 -15 -10 -5 0 5 10 15 20 



STERN PUNE ANGLE 

Figure 17. Steady State Roll Angle 



PITCH angle theta vs DS 




STERN PUNE angle 

Figure 18. Steady State Pitch Angle 
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It can be seen that all horizontal plane variables, v, p, r and 4», exhibit the typical 
characteristics of pitchfork bifurcation. For 6r= 0, there exists a dive plane angle 6s 
between -5 and -10 degrees, such that the in plane solution becomes unstable and 
symmetric stable out of plane solutions appear. Nonzero values of 5r result, naturally, 
in out of plane solutions for the entire range of 6s. Multiple solutions appear also at a 
certain value of 6s, called turning point, which is a function of 6r. 

The vertical plane variables u, w, q, and 6 are even functions of horizontal plane 
motions and, as a result, they exhibit "solution separation" of the in plane motions for 
values of 6s lower than the bifurcation point. Of special interest is an examination of the 
steady state values of the pitch angle 6 and roll angle since these are directly related 
to the vehicle orientation. A clear examination of the possible steady state solution, so 
is the pair (v-6, t-#). When 6r= 0, it can be seen from Figure 18 that the stable in- 
plane solution exhibits a steady state pitch angle greater than 90 degrees for values of 6s 
less than about -5 degrees. This is, of course, true up to the pitchfork bifurcation point. 
This phenomenon of inverted pendulum stabilization was first discussed by McKinley 
[Ref. 14]. It appears, however, from the results of Figure 18 that the inverted pendulum 
stabilization is highly degenerate and is destroyed as soon as a non-zero rudder angle 6r 
is applied. The two solutions, 6 and r-9, do not cross each other as they did for 6r=0. 
Instead, they veer of each other, the lower one corresponding to ^=0 and the upper one 
to 4 >=t. Asa result, the vehicle configuration is always keel down, the only difference 
between the two solution sets at Figures 17 and 18 being 180 degrees heading. 
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Another important observation can be made by considering the results of Figures 
16 and 17, in terms of the roll angle It appears that there exists a range of stemplane 
angle where the exact value of is very sensitive to the exact value of 6r for small 
values of 6r, The rudder angle is related to the angle of attack of the fluid flow, which 
is highly dq)endent on such factors as vortex shedding, flow separation, and ambient 
flow turbulence Saipkaya [Ref. 5]. Since these factors are to a great extent uncertain and 
difficult to model, situations like the one described above should be avoided in practice. 
This bifurcation and continuation methods used in this thesis provide a systematic and 
effective way of identifying precisely parameter regions where such sensitivity of the 
response to parameter values exists. These regions should be avoided when executing a 
proper recovery procedure. 
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rv. CONTINUATION IN 8r 



A. GENERAL 

In this chapter, the results of continuation in terms of the rudder angle will be 
discussed. In this continuation run, the initial conditions were taken from the previous 
results of continuation in terms of the rudder angle. The configuration of the vehicle was 
just like the previous case. In this case the system of equations (35) was solved for 
changing values of the primary parameter X (here X= 6r). At the end of the runs, a 
hysteresis behavior was observed for some of the state variables. This hysteresis 
behavior, which is complementary to pitchfork bifurcation, is explained in the following 
section. 

B. HYSTERESIS PHENOMENON 

When one-parameter families of equation are studied, there is no need to talk about 
a hysteresis phenomenon. But similarly to the pitchfork bifurcation, this situation might 
change when a two-parameter family of equations is studied, 

f(y,X,a) . 

Now let us consider the following equation: 
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x^-A, + ax = 0 



(41) 



Let us call this algebraic equation as a function of x F(x) = 0. It can be observed that 
equation (41) is similar to the generic pitchfork equation (40). The only difference is that 
the roles of X and a have been interchanged. The constant term is now our primary 
bifurcation parameter, while the coefficient of x is the perturbation parameter. We know 
that at the turning point F, and dF/dx, should be zero. So, 



dx 



=0 



x^ = - 



oc 

3 



If we substitute this value into Equation (41), we will get the following equation. 

27A,2+4a3=0 (42) 



This equation gives the relationship between X and a at the vertical tangent. Now we 




Figure 19. Cusp Curve 



39 



consider equation (41) as a two parameter model with X and a having equal rights. 
Equation (42) establishes a curve in the (a , X) plane, which takes the form of a cusp 
(Figure 19), (Cusp curves will be explained in a more general context in Chapter Vin.). 

For all combinations of ( a, X) values inside the cusp (hatched region), Equation 
(41) has three solutions; for values outside the cusp, there is only one solution. This 
becomes clear when branching diagrams are drawn. As a generic behavior, if we draw 
branching diagrams for different values of a when X is the primary continuation 
parameter, we get the Figures 20, 21, 22 according to the sign of a. 

• a = 0 

• a > 0 

• a < 0 

X 




Figure 20. Branching Diagram When a = 0 
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Figure 21. Branching Diagram When a < 0 




In this last case, we see a hysteresis effect. When hysteresis occurs in a branch, 
there are two associated turning points with hysteresis behavior. These points occur 

I 

precisely at the two real solutions of X from Equation 42. Here there is a jump when the 
solution is on an appropriate part of the branch. The two outer solution in x are stable 
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while the inner most solution between the two turning points is unstable. The two dotted 
straight lines of Figure 19 represent two distinct paths through the cusp and explained 
in Chapter Vm. 

C. STEADY STATE RESULTS 

Figures 23 through 30 show steady state solutions for the eight state variables, 
namely; u, v ,w, p, q, r, B. In Figures 23 through 30 solid lines correspond to 6s= 
0 degree, dashed lines to 5s = -5 degrees, dotted lines to 6s = -10 degrees and dash-dot 
lines to 6s = -20 degrees. For these continuation runs, the rudder angle (6r) is the 
primary continuation parameter (i.e,X) and stem plane angle (6s) is perturbation 
parameter (i.e,a) and values of all other parameters are kept fixed:excess buoyancy, 
6B = 2 % of the vehicle weight (W); deflection of bow planes, 6b = 0 ; location of 
horizontal and vertical centers of buoyancy ,Xb= Zg = 0;location of longitudinal center 
of buoyancy, x^g =-l % of the vehicle length (L). 
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SURGE VELOCITY U vs DR 



z> 




Figure 23. 



Steady State 



RUDDER angle 
Surge Velocity 



SWAY VELOCITY V vs DR 




Figure 24. 



Steady State 



RUDDER angle 
Sway Velocity 
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HEAVE velocity W vs or 




ROLL ANGULAR VELOCITY P vs DR 




RUDDER angle 

Figure 26. Steady State Roll Angular Velocity 
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0.03 



PITCH angular velocity U vs dr 




YAW angular velocity R vs dr 




Figure 28. 



rudder angle 

Steady State Yaw Angular Velocity 
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THETA 



ROLL angle phi vs dr 




Figure 29. 
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Steady State Roll Angle 
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Figure 30. 



RUDDER ANGLE 

Steady State Pitch Angle 



46 




Figures 23 through 30 show that the horizontal plane state variables v, p, r, 
behave like Figure 22 and exhibit a hysteresis behavior for certain values of 6s. If we 
zoom in Figure 28, we obtain Figure 31 and can see this hysteresis effect very clearly. 
When 6s = -15 degrees, there is a typical hysteresis phenomenon. Figures 24, 26, and 
28 also show the same kind of behavior. 



TYPICAL hysteresis behavior from continuation runs 




RUDDER angle 

Figure 31. Typical Hysteresis Effect from Continuation 

On the other hand, the vertical plane state variables u, w, q, 6 behave like an even 
function in x (for example | x | ). This behavior can be seen in Figures 23, 25, 27, 30. 
where it is clear that they exhibit generic behavior of an even function for different 
values of the pertuibation parameter (i.e, 6s). If we zoom in Figure 23, we can see this 
behavior, very clearly. 
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TYPICAL BEHAVIOR OF AN EVEN FUCTION FROM CONTINUATION 




RUDDER angle 

Figure 32. An Even Function from Continuation 
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V. CONTINUATION IN x, 



GB 



A. GENERAL 

In this chapter, the results of continuation in terms of the longitudinal center of 
buoyancy will be discussed. In continuation runs, the system of equations (35) was 
solved for changing values of the primary parameter X where in this case X = x<-b). In 
these runs, first we fixed the value of rudder angle and changed the values of stem plane 
angle, then we kept the stemplane angle fixed and changed the rudder angle, to see how 
the steady state values were changing. When we keep the rudder angle at zero and 
change the values of stemplane angle, we observe that there is a pitchfork bifurcation for 
the states v, p, r. Then, we fixed 6s at -20 degrees and started to change mdder angle 
to see how the pitchfork bifurcation diagrams were changing with changing mdder angle. 
As a last continuation mn for this configuration, the value of stemplane angle was fixed 
at zero and we changed the value of mdder angle. We didn’t see any solution branching 
for this case. 

B. STEADY STATE RESULTS 

1. 6r=0 Degrees, 6s is Perturbation Parameter 

Figures 33 through 40 show steady state solutions for eight state variables, 
namely; u, v, w, p, q, r, $. In Figures 33 through 40, solid lines correspond to 6s=0 
degrees, dotted lines correspond to 6s=-5 degrees, and dash-dot lines correspond to 6s=- 
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20 degrees. In graphs of the state, variables v, p, r, d dash-dot lines may seem like solid 
lines; this is because the continuation between the parameter values Xqb = -13 degrees 
and Xqb = .7 degrees traced the same curve twice in order to make sure that all solutions 
for all the state variables were captured. Therefore, out of plane solutions in this 
parameter range belong to stemplane angle 6s = -20 degrees. 

In all these cases, the longitudinal center of buoyancy is the continuation parameter 
and 6s is the perturbation parameter, while 6r is kept at zero degrees. Values of the other 
parameters are kept fixed: excess buoyancy, 6B =2 % of the vehicle weight (W); 
deflection of the bow planes, 6b = 0: location of horizontal and vertical centers of 
gravity, Xg = Zb =0: location of vertical center of gravity, Zgb= -1 feet. 



SURGE velocity U vs XGB 




LONGITUDINAL CENTER OF BUOYANCY %L 

Figure 33 . Steady State Surge Velocity 
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SWAY VELOCITY V vs XGB 
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Figure 34. Steady State Sway Velocity 



HEAVE VELOCITY W vs XGB 




longitudinal center of buoyancy %l 

Figure 35. Steady State Heave Velocity 
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ROLL ANGULAR VELOCITY P vs XGB 
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Figure 36* Steady State Roll Angular Velocity 
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Figure 37. Steady State Pitch Angular Velocity 



PITCH ANGULAR VELOCITY Q vs XGB 
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YAW ANGULAR VELOCITY R vs XGB 
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Figure 38. Steady State Yaw Angular Velocity 
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Figure 39. 



ROLL ANGULAR VELOCITY PHI vs XGB 
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PITCH ANGUUVR VELOCITY THETA vs XGB 




2. 5s=-20 Degrees, ir is the Perturbation Parameter 

Figures 41 through 48 show steady state solutions for the state variables u, 
V, w, p, q, r, 4>, 6 . In Figures 41 through 48, solid lines correspond to 5r=0 degrees, 
dotted lines correspond to 6r=2.5 degrees, and dash-dot lines correspond to 6r=5 
degrees. 

In all of these cases, longitudinal center of buoyancy is the primary 
continuation parameter and 6r is the perturbation parameter, while 6s is kept at -20 
degrees. Values of the other parameters are kept fixed: Excess buoyancy, 6B = 2 % of 
the vehicle weight (W): deflection of bow planes, 6b = 0; location of horizontal and 
vertical centers of buoyancy, Xg = Zg =0; location of the vertical center of gravity, Zo= 
.1 feet. 
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SURGE velocity U vs XGB 




longitudinal center of buoyancy %l 

Figure 41. Steady State Surge Velocity 



SURGE velocity V vs XGB 




LONGITUDINAL CENTER OF BUOYANCY ZL 

Figure 42. Steady State Sway Velocity 
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HEAVE VELOCITY W vs XGB 




ROLL ANGULAR VELOCITY P vs XGB 




longitudinal CENTER OF BUOYANCY JSL 

Figure 44. Steady State Roll Angular Velocity 
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PITCH angular velocity Q vs xgb 




LONGITUDINAL CENTER OF BUOYANCY JSL 

Figure 45. Steady State Yaw Angular Velocity 



YAW ANGULAR VELOCITY R vs XGB 




LONGITUDINAL CENTER OF BUOYANCY ZL 

Figure 46. Steady State Yaw Angular Velocity 
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THETA 



YAW angle phi vs XGB 




Figure 47 . Steady State Roll Angle 



PITCH angle theta vs XGB 




LONGITUDINAL CENTER OF BUOYANCY 55L 

Figure 48. Steady State Pitch Angle 
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3. 5s = 0 Degrees, 6r is the Perturbation Parameter 

Figures 49 through 56 show steady state solutions for the eight state 
variables, namely u, v, w, p, q, r, 4>, 6 . In Figures 49 through 56, solid lines 
correspond to6r = 0 degrees, dashed lines correspond to 6r = 0.5 degrees, dotted lines 
correspond to 6r = 2.5 degrees and dash-dot lines correspond to 5r = 5 degrees. 

In all these cases, the longitudinal center of buoyancy is the primary 
continuation parameter and 6r is the perturbation parameter. Values of the other 
parameters are kept fixed: Stemplane action, 6s = 0 degrees: excess buoyancy, SB = 2 
% of the vehicle weight (W); deflection of bow planes, 6b = 0; location of horizontal 
and vertical centers of buoyancy, Xg = Zg = 0; location of vertical center buoyancy, Zqb 
= 0.1 feet. 



SURGE VELOCITY u vs XGB 




longitudinal center of buoyancy %L 
Figure 49 . Steady State Surge Velocity 
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SWAY VELOCITY V vs XGB 




HEAVE velocity W vs XGB 




longitudinal center of buoyancy %l 

Figure 51. Steady State Heave Velocity 
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ROLL ANGULAR VELOCITY P vs XGB 




longitudinal center op buoyancy ZL 

Figure 52. Steady State Roll Angular Velocity 



x10-3 PITCH ANGULAR VELOCITY Q vs XGB 




LONGITUDINAL CENTER OF BUOYANCY ZL 

Figure 53. Steady State Pitch Angular Velocity 
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YAW ANGULAR VELOCITY R vs XGB 
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Figure 55. 
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THETA 



PITCH angle theta vs XGB 




LONGITUDINAL CENTER OF BUOYANCY ZL 

Figure 56. Steady State Pitch Angle 
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VI. CONTINUATION IN 5B 



A. GENERAL 

In this chapter, the results of continuation in terms of the excess buoyancy will be 
discussed. In continuation runs, the system of equations (35) was solved for changing 
values of the primary parameter X (in this case X= 6B). In these runs, first we fixed the 
value of rudder angle and changed the values of stemplane angle, then we kept the 
stemplane angle fixed at zero and changed the rudder angle to see how the steady state 
values were changing. When we kept the rudder angle at zero and changed the values of 
stemplane angle, we saw that there was a pitchfork bifurcation for the state variables v, 
p, r, and ^». After we observed this pitchfork bifurcation for 5s= -20 degrees and 6r=0 
degrees, we started to change the rudder angle 6r to see how the pitchfork bifurcation 
diagrams were changing. As an other continuation mn we fixed 6r=0 degrees and 
changed 6s towards positive values; however, we didn’t see any branching in this case. 
As a last continuation mn for this configuration, the value of stemplane angle was fixed 
at zero and the value of mdder angle was changed. We didn’t see any branching behavior 
for this case. 
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B. STEADY STATE DIAGRAMS 



1. ir=0 Degrees, 5s is the Perturbation Parameter 

Figures 57 through 64 show steady state solutions for eight state variables, 
namely; u, v, w, p, q, r, 4», 6 . In Figures 57 through 64, solid lines correspond to 6s =0 
degrees, dashed lines correspond to 6s=-5 degrees, and dash-dot lines correspond to 
6s =-20 degrees. In all of these cases, excess buoyancy is the continuation parameter, 
while 6r is kept at zero degrees. Values of the other parameters are kept fixed: 
longitudinal center of buoyancy Xcb=- 1 % of the vehicles length (L); deflection of the 
bow planes, 6b =0; location of horizontal and vertical centers of buoyancy, Xb=Zb= 0; 
location of vertical center of buoyancy, Zgb= 1 feet. 

As we can see from Figures 57 through 64, there is a pitchfork bifurcation for state 
variables u, v, p, r, 4>. We have out of plane solutions for stemplane angle 6s =-20 
degrees, this is also the configuration in which we observe bifurcations. 
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SURGE velocity U vs DELS 




EXCESS BUOYANCY J5W 
Figure 57. Steady State Surge Velocity 
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SWAY velocity V vs DELB 




0 0.5 1 1.5 2 2.5 3 



EXCESS BUOYANCY %W 

Figure 58. Steady State Sway Velocity 



HEAVE VELOCITY W vs DELB 




EXCESS BUOYANCY %W 

Figure 59. Steady State Heave Velocity 
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ROLL angular velocity P vs DELB 
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Figure 60. Steady State Roll Angular Velocity 
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YAW angular velocity R vs DELB 
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Figure 62. 



EXCESS BUOYANCY %W 

Steady State Yaw Angular Velocity 
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ROLL angle phi vs delb 




EXCESS BUOYANCY ZW 

Figure 63. Steady State Roll Angle 
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PITCH angle theta vs delb 




After analyzing the above diagrams, we decided to look at a case with a 
positive stemplane angle. We observed that there was no solution branching for this 
case. We got the inplane solutions and there were no out of plane solutions. In Figure 
65, we can see the situation for the state variable u. In this figure, solid lines correspond 
to 6s=0 degrees, dotted lines correspond to 6s =5 degrees, and dash-dot lines correspond 
to 6s = 20 degrees. The steady state solutions are not very sensitive to positive stemplane 
action. 
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0 0.5 1 1.5 2 2.5 2 

EXCESS BUOYANCY J5W 

Figure 65. Steady State Surge Velocity When 6s is Positive 

2. 6s = -20 Degrees, 5r is the Perturbation Parameter 

Figures 66 through 73 show steady state solutions for the state variables u, 
V, w, p, q, r, 6 . In these figures, dashed lines correspond to 6r= 0.5 degrees, and 
dash-dot lines correspond to 6r= 5 degrees. In all of these the cases, excess buoyancy 
is the primary continuation parameter and 5r is the perturbation parameter, while 6s is 
kept -20 degrees. Values of the other parameters are kept fixed. Longitudinal center of 
buoyancy, -1 % of the vehicle length (L); deflection of bow planes, 6b= 0; location of 
horizontal and vertical centers of buoyancy, Xb= z^= 0; location of the vertical center 
of buoyancy, Zqb= 0.1 feet. By doing this, we saw how the bifurcation points were 
changing as we changed the perturbation parameter 6r. 
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SURGE VELOCITT U vs DELB 




EXCESS BUOYANCY %'N 

Figure 66. Steady State Surge Velocity 



SWAY velocity V vs DELB 




Figure 67 . 



Steady State 



EXCESS BUOYANCY 7,'N 

Sway Velocity 
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HEAVE VELOCITT W vs DELB 




EXCESS BUOYANCY nw 

Figure 68. Steady State Heave Velocity 



ROLL ANGULAR VELOCITY P vs DELS 




EXCESS BUOYANCY J5W 

Figure 69 . Steady State Roll Angular Velocity 
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0.035 



PITCH angular velocity Q vs DELB 




EXCESS BUOYANCY JSW 

Figure 70. Steady State Pitch Angular Velocity 



YAW ANGULAR VELOCITY R vs DELB 




EXCESS BUOYANCY J5W 

Figure 71. Steady State Yaw Angular Velocity 
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THETA 



ROLL angle phi vs OELB 




EXCESS BUOYANCY J5W 

Figure 72. Steady State Roll Angle 



PITCH angle theta vs delb 




EXCESS BUOYANCY %W 

Figure 73. Steady State Pitch Angle 
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3. 5s =0 Degrees, 5r is Perturbation Parameter 

Figures 74 through 81 show steady state solutions for the eight state 
variables, namely u, v, w, p, q, r, 4», 6. In Figures 74 through 81 solid lines correspond 
to 6r= 0 degrees, dashed lines correspond to 6r= 0.5 degrees, and dash-dot lines 
correspond to 6r=5 degrees. In aU these cases, the excess buoyancy 6B is the primary 
continuation parameter and 6r is the perturbation parameter. Values of the other 
parameters are kept fixed: Stemplane action, 6s =0 degrees; longitudinal center of 
buoyancy Xgb=- 1% of the vehicle length (L); deflection of bowplanes, 6b =0; location 
of horizontal and vertical centers of buoyancy, Xb= Zb= 0; location of vertical center of 
buoyancy, Zqb— 01 As we can see, there was no solution branching for this case. 
Of course, we got out of plane solutions for nonzero rudder angle values. 



SURGE velocity U vs DELB 




EXCESS BUOYANCY J!W 

Figure 74. Steady State Surge Velocity 
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SWAY VELOCITY V vs DELB 




HEAVE VELOCITY W vs DELB 




EXCESS BUOYANCY nw 

Figure 76. Steady State Heave Velocity 
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0.025 



ROLL ANGULAR VELOCITY P vs DELB 




x10-3 PITCH angular velocity Q vs DELB 




EXCESS BUOYANCY 55W 

Figure 78. Steady State Roll Angle 
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0.03 



YAW ANGULAR VELOCITY R vs DELB 




Figure 79. 
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Steady State Yaw Angular Velocity 
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ROLL angle phi vs delb 
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Figure 80. Steady State Roll Angle 
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THETA 



PITCH angle theta vs DELB 




EXCESS BUOYANCY %VJ 

Figure 81. Steady State Pitch Angle 
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Vn. CONTINUATION IN z, 



GB 



A. GENERAL 

In this chapter, the results of continuation in terms of the vertical center of 
buoyancy will be discussed. In continuation runs, the system of equations (35) was solved 
for changing values of the primary parameter X (in this case X= Zgb)- In these runs, first 
we fixed the value of bidder angle and changed the values of stemplane angle, then we 
kept the stemplane angle fixed at zero and changed the mdder angle to see how the 
steady state values were changing. When we kept the rudder angle at zero and changed 
the values of stemplane angle, we saw that there was a pitchfork bifurcation for the state 
variables u, v, w, p, r, 6 . After we saw that there was a pitchfork bifurcation when 
5s = -20 degrees and 6r= 0 degrees, we started to change mdder angle 6r to see how the 
pitchfork bifurcation diagrams were changing. As a last continuation mn for this 
configuration, the value of stem plane angle was fixed at zero and the value of mdder 
angle was changed. We didn’t observe any branching behavior for this case. 

B. STEADY STATE RESULTS 

1. 5r= 0 Degrees, 5s is the Perturbation Parameter 

Figures 82 through 89 show steady state solutions for eight state variables, 
namely; u, v, w, p, q, r, 6 . In Figures 82 through 89, solid lines correspond to 5s = 
0 degrees, dashed lines correspond to 5s = -5 degrees and dash-dot lines correspond to 
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5s =-20 degrees. In all of these cases, the vertical center of buoyancy is the primary 
continuation parameter, and 5s is the perturbation parameter, while 5r is kept at zero 
degrees. Values of the other parameters are kept fixed: Excess buoyancy 5B=2% of the 
vehicle weight (W); longitudinal center of buoyancy Xgb=- 1% of vehicle length; 
deflection of the bow planes 5b =0; location of the vertical and horizontal centers of 
buoyancy, Zb=Xb= 0. As we see from Figures 82 through 89, there is a pitchfork 
bifurcation for state variables u, v, w, p, r, 4>, 0 . We have out of plane solutions for 
this vehicle configuration. 



SURGE VELOCITY U vs ZGB 




VERTICAL CENTER OF BUOYANCY (f«et) 
Figure 82. Steady State Surge Velocity 
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SWAY velocity V vs ZGB 
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Figure 83 . Steady State Sway Velocity 



HEAVE VELOCITY W vs ZGB 
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Figure 84. Steady State Heave Velocity 
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ROLL ANGULAR VELOCITY P vs ZGB 
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Figure 85. Steady State Roll Angular Velocity 



0.3 



0.016 
0.014 
0.012 
0.01 
O 0.008 
0.006 
0.004 
0.002 
0 



PITCH ANGULAR VELOCITY Q vs ZGB 





1 .J 

✓ 

/ 


* ^ 

* s 










/ 

/ 

/ 




/ 

/ 

/ 








J 




S 

\ 

.'n 






/ 

t 

1 








\ 

\ 

\ 




/ 

/ 

* 








\ 

\ 

\ 




f 

/ 








\ 

\ 




f 










\ 

\ 

\ 












V 

\ 

\ 



0.05 



0.1 



0.15 



0.2 



0.25 



0.3 



Figure 86. 
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Steady State Pitch Angular Velocity 
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YAW ANGULAR VELOCITY R vs ZGB 
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Figure 87. Steady State Pitch Angular Velocity 
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Figure 88. 
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PITCH angle theta vs ZGB 




VERTICAL CENTER OF BUOYANCY (fe«t) 

Figure 89. Steady State Pitch Angle 

2. 5s =-20 Degrees, 5r is the Perturbation Parameter 

Figures 90 through 97 show steady state solutions for the state variables u, 
V, w, p, q, r, 6 . In Figures 90 through 97, solid lines correspond to 6r=0 degrees, 
dashed lines when 5r=0.5 degrees, dash-dot lines correspond to 6r=5 degrees. In all 
of these cases, vertical center of buoyancy Zqb is primary continuation parameter and 5r 
is the perturbation parameter, while 6s is kept at -20 degrees. Values of the other 
parameters are kept fixed; excess buoyancy 6 b= 2% vehicle weight (W), longitudinal 
center of buoyancy x^b— - 1 % of the vehicle length (L); deflection of bowplanes 6b =0; 
location of horizontal and vertical centers of buoyancy, Xb=Zb= 0. By doing this, we 
saw how the bifurcation points were changing as we change perturbation parameter 6r. 
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A bifurcation phenomenon took place for 5r=0 degrees and disappeared after we applied 
a small nonzero rudder angle value. 



SURGE VELOCITT U vs ZGB 




VERTICAL CENTER OF BUOYANCY (feet) 

Figure 90. Steady State Surge Velocity 



SWAY VELOCITY V vs ZGB 
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Steady State Sway Velocity 
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Figure 91 



HEAVE VELOCITY W vs ZGB 




ROLL ANGULAR VELOCITY P vs ZGB 




vertical center or buoyancy (feet) 

Figure 93 . Steady State Roll Angular Velocity 
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0.025 



PITCH ANGULAR VELOCITY Q vs 2GB 




VERTICAL CENTER OF BUOYANCY (feet) 

Figure 94. Steady State Pitch Angular Velocity 



YAW ANGULAR VELOCITY R vs ZGB 
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Figure 95. Steady State Yaw Angular Velocity 
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THETA 



ROLL angle phi vs ZGB 




PITCH angle theta vs ZGB 




vertical center of buoyancy (fe«t) 

Figure 97 . Steady State Pitch Angle 
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3. 5s = 0 Degrees, 5r is the Perturbation Parameter 

Figures 98 through 105 show steady state solutions for the eight state 
variables, namely u, v, w, p, q, r, B. In Figures 98 through 105, solid lines 
correspond to 6r=0 degrees, dashed lines correspond to 6r=.5 degrees, dash-dot lines 
correspond to 6r=5 degrees. In all these cases, the vertical center of buoyancy is the 
primary continuation parameter. Values of the other parameters are kept fixed; stemplane 
action 6s =0 degrees; longitudinal center of buoyancy Xqb = -1 % of vehicle length (L); 
deflection of mow planes, 5b =0 degrees; excess buoyancy 6B= 2% of vehicle weight 
(W); location of horizontal and vertical centers of buoyancy Xb=Zb=0. As we see, there 



SURGE VELOCITY U vs ZGB 




vertical center of buoyancy (feet) 
Figure 98. Steady State Surge Velocity 
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is no solution branching for this last case. We have out of plane solutions for nonzero 
rudder angles. 



SWAY velocity V v» ZG0 




VERTICAL CENTER OF BUOYANCY (f««t) 

Figure 99. Steady State Sway Velocity 
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Figure 100. Steady State Heave Velocity 



HEAVE VELOCITY W vs ZGB 
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Figure 101. Steady State Roll Angular Velocity 



x10-3 PITCH ANGULAR VELOCITY Q vs ZGB 
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Figure 102. Steady State Pitch Angular Velocity 
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YAW angular velocity R vs zgb 



0 

- 0.005 

- 0.01 

- 0.015 

- 0.02 

- 0.025 

- 0.03 

- 0.035 

- 0.04 



Figure 103 . 
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Figure 104. Steady State Roll Angle 
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THETA 



PITCH angle theta vs ZGB 
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Figure 105. Steady State Pitch Angle 
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vm. CUSP CURVES 



A. GENERAL 

In this chapter, the cusp phenomenon will be discussed. After seeing the branching 
characteristics of the steady state solutions of our system, we decided to look at how the 
bifurcation points and the turning points were changing with changing values of rudder 
angle. The main difference between a bifurcation point and a turning point is that, at a 
bifurcation point, exactly two branches intersect with two distinct tangents; whereas, at 
a turning point, the branch comes from one side and turns back. We can see this 
phenomenon from Figure 106. This is the graph that was generated for the surge velocity 
when the primary continuation parameter was XcB^nd the stemplane angle 6s =-20 
degrees. In Figure 106, points A and B are bifurcation points, and points C and D are 
turning points. Clearly, locally there are no solutions on one side of a turning point and 
two solutions on the other side, but there are solutions on each side of a bifurcation 
point. 

Let us look at Equation 43 again, this is the model equation that represents a cusp, 

x^-Bx+A=Q ( 43 ) 

Let us call this algebraic equation as a function of x, F(x)=0. For the critical curve, F 
and dF/dx should be zero. 
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TURNING ANO BIFURCATION POINTS 




Figure 106. Turning and Bifurcation Points 
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Figure 107 . Cusp Curve 
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If we substitute this value into equation (43), we get the following relation; 

If we plot this relation in the AB plane, we get the cusp curve shown in Figure 107. This 
represents a sq)aration. Equation 43 has three solutions for parameter values inside the 
cusp and one solution outside this parameter values. In Figure 107, if we fix the value 
of A and vary the value of B (path (a) in Figure 107) we get a pitchfork bifurcation 
diagram, as shown in Figure 108. 



constant a, changing b 




Figure 108. Pitchfork Bifurcation 

In Figure 107, if we fix the value of the value of B and vary the value of A (path (b) in 
Figure 107) we get a hysteresis behavior, as shown in Figure 110. It is clear, that when 
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a line crosses the cusp curve in the (A, B) plane a branching phenomenon takes plane. 
As mentioned before, Equation 43 has three solutions for parameter values inside the 
cusp and one solution outside the cusp this can be seen easily in Figure 109. Now let us 
take a look at Figure 111. As can be seen, there is a simple way of recovering a 
bifurcation diagram from the path through the given cusp by lifting the path to the cusp 
surface. 




Figure 109. 3-D Graph of Equation 43 
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Figure 110. Hysteresis Diagram 




Figure 111. Geometry of Cusp Phenomenon 
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If we have the cusp curves for a multiparameter equations system, we can always 
have an idea about the solutions of this system for different perturbed cases. This means 
that by looking at the path of the parameters through the cusp curve, we can see how 
many solutions we have for the corresponding configuration of the system. Similarly, we 
can think about another problem which has a turning point. If we draw the position of 
turning point in the plane with changing parameter values, we will get two solutions 
inside the curve and no solutions outside these parameter values. Therefore, these graphs 
can tell us the nature of the solution set of our system. In our study, we tried to get the 
cusp curves and the position of the turning points and the results are presented in the 
following section. 

B. RESULTS 

First of all, let us look at the case when we have the stemplane action as our 
primary continuation parameter. When 6s is between -5 and -10 degrees, we had a 
bifurcation point and in that case the rudder angle was zero, see Figure 13. To get the 
complete cusp curve, we changed the value of 6r and found the locations of the 
bifurcation point. These were plotted as 6s versus 6r, see Figure 112. 
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CUSP CURVE WHEN CONTINUATION PARAMETER is DS 




Figure 112 . Cusp Curve When Continuation Parameter is 

Sternplane angle 



In Chapter V, we saw that when we have Xcb as primary continuation parameter, 
and sternplane action 6s=-20 degrees, we have both bifurcation and turning points. 
Figure 106. To see the effect of changing the rudder angle on the bifurcation and turning 
points, we changed the value of rudder angle and found the location of bifurcation and 
turning points in the (x^b, dr) plane. The results are shown on the following graphs. The 
bifurcation points generated a cusp curve and the turning points generated a parabola in 
the (xBo.dr) plane. 
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CUSP CURVE 




XGB 

Figure 113 . Cusp Curve When Continuation Parameter is x^g 



TURNING POINT 




XGB 

Figure 114. Turning Point When Continuation Parameter is x^g 
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CUSP CURVE 







turning point 




X08 

Figure 116. Second Turning Point When Continuation Parameter 
is Xqb 
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In Chapter VI, we used the excess buoyancy as our primary continuation 
parameter. In that configuration, we had a turning point and a bifurcation point for 
stemplane angle 6s =-20 degrees, (Figure 57 in Chapter VI). To see the effect of 
different rudder angles on these points, we changed the value of the rudder angle and 
found the location of these points in the (Xqb, 6r) plane. These results are shown in 
Figure 117 and Figure 118. Again, we got a cusp curve and a parabola. 



CUSP CURVE 




delta B 

Figure 117. Cusp Curve When Continuation Parameter Is Excess 
Buoyancy 
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TURNING POINT 




delta B 

Figure 118. Turning Point When Continuation Parameter is 
Excess Buoyancy 



In Chapter Vn, we used the vertical center of buoyancy as our primary 
continuation parameter. In that continuation, run we found that there were two 
bifurcation points at stemplane angle 6s =-20 degrees, see Figure 82 Chapter Vn. In that 
figure, the bifurcation point when Zqb is around zero is not clear, but it becomes clear 
when we changed the value of rudder angle. To construct the cusp curves, we changed 
the value of rudder angle and found the coordinates of bifurcation points in the (Zqb, 6r) 
plane. These graphs are shown in Figures 119 and 120. 
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CUSP CURVE 




ZGB 

Figure 119. Cusp Curve When Continuation Parameter is Zqq 



CUSP CURVE 




Figure 120. Second Cusp Curve When Continuation Parameter is 

^GB 
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As a final study in this chapter, we changed the value of ya in the first continuation 
run and tried to get the biased cusp curves for this perturbed case. As a result, we can 
still see a symmetric bifurcation for this perturbed case for a specific parameter 
combination. It is clear, however, that this situation is very degenerate and can easily be 
destroyed for different parameter values. In Figure 121, solid lines correspond to yc=0 
feet, dashed lines correspond to yc=0.001 feet, dotted lines correspond to yc =0.003 
feet, and dash-dot lines correspond to yo =0.005 feet. 




DS 

Figure 121. Cusp Curve When Continuation Parameter is 6s for 
Different Values of 
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K. DYNAMIC STABILITY ANALYSIS 



A. GENERAL 

So far, we have performed continuation runs and have obtained steady state 
solutions of the system. For stability analysis, we chose a representative case and 
performed the stability analysis for that case. We selected the case when we had Xqb as 
our primary continuation parameter and rudder angle 6r=0 degrees, and stemplane angle 
6s=-20 degrees. To predict dynamic stability, we used the six degree of freedom 
equations of motion along with the Euler angle rate equations for the derivatives of the 
angles of pitch and roll (0 and 4 ). These equations of motion were then linearized around 
the steady state nominal points computed in the previous chapters. Eigenvalue analysis 
was then used to compute the stability characteristics of each solution. 

B. DYNAMIC STABILITY RESULTS 

The stability analysis results are shown in Figure 122 and Figure 123. In these 
figures, solid lines correspond to stable solutions while dashed lines correspond to 
unstable solutions. In Figure 122 up to point A, we have a stable solution, which is 
predominantly forward motion, because we have a high surge velocity and the vehicle 
is moving in the vertical plane. Between point A and B an out of plane solution exists 
and the vehicle exhibits motion in the horizontal plane, also. From Figure 123, it is clear 
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that at point A the sway velocity becomes nonzero. This out of plane solution is stable 
up to point B, where 



SURGE VELOCITV 




Figure 122 . Steady State Surge Velocity 



SWAY VELOCITY 




Figure 123. Steady State Sway Velocity 
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the out of plane motion becomes unstable up to point C, where it disappears. In Figure 
122, the motion of the vehicle from point D to E is stable and has a low surge velocity, 
which means that the motion of the vehicle is predominantly in the upward direction. 
The vehicle is moving in the vertical plane and there is no horizontal plane motion. 

To assess stability of each solution, we used eigenvalue analysis. The eigenvalues 
that we computed for the following cases are summarized in Table 1 . 

• Case 1 

Stemplane angle 6s=-20 degrees, rudder angle 6r=0 degrees, Xqb=- 2 % L, in 
plane solution. 

• Case 2 

Stemplane angle 6s=-20 degrees, rudder angle 6r=0 degrees, Xgb=-1.25 % L, out 
of plane solution. 

• Case 3 

Stemplane angle 6s=-20 degrees, rudder angle 6r=0 degrees, Xob= 0.57% L point 
C Figure 122. 

• Case 4 

Stemplane angle 6s=-20 degrees, mdder angle 6r=0 degrees, x<3b=0 % L. This 
is oh the stable branch between point D and £ in Figure 122. 

Table 1 represents the eigenvalues of the system. We can see for case three that 
there is a pair of complex conjugate eigenvalues with positive real parts. This means that 
we have a hopf bifurcation prior to that point, and we expect periodic solutions for that 
part of the solution branch. Since we have a stable solution up to point E in Figure 122, 
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we don’t expect to see this periodic out of plane behavior up to point E. The methods 
that we used in this work give only stationary solutions, therefore, to see this periodic 
behavior we used numerical simulation techniques. 

TABLE 1 



Eigenvalues 


Case 1 


Case 2 


Case 3 


Case 4 




-0.962-1-0.499 


-2.407 


-2.591 


-0.093-1-0.522 


K 


-0.962-0.499 


-0.855-1-0.413 


-0.968 -»- 0.427 


-0.093-0.522 


X3 


-0.068-1-0.014 


-0.855-0.413 


-0.968-0.427 


-0.067-1-0.064 


X4 


-0.068-0.014 


-0.723 


-1.048 


-0.067-0.064 


X3 


-0.761 


-0.066 


-0.314 


-0.193-1-1.531 


K 


-2.731 


-0.018-1-0.044 


-0.129 


-0.193-1.531 


X7 


-0.198-h0.115 


-0.181-0.044 


-1-0.015-1-0.004 


-0.240 


^8 


-0.198-0.115 


-0.271 


-fO.015-0.004 


-0.045 



C. NUMERICAL INTEGRATION RESULTS 

The linearized dynamic response results were verified by simulations using 
numerical integration of the full six degrees of freedom equations of motion for the 
swimmer delivery vehicle (SDV). Figures 124 and 125 show a plot of surge velocity 
(U), and sway velocity (V) versus time respectively for the center of gravity aft of the 
center of buoyancy case (Xgb=-2) with a diveplane angle (6s) of -20 degrees, and rudder 
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SURGE velocity 




Figure 124. Time History of Surge Velocity for Case 1. 



SWAY VELOCITY 




Figure 125. Time History of Sway Velocity for Case 1. 
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angle (6r) of 0 degrees. The steady state results of the numerical integration method 
match the linearized dynamic results exactly. 

Figures 126 and 127 show a plot of surge velocity (U), and sway velocity (V) 
versus time respectively for the center of gravity aft of the center of buoyancy case 
(xgb= -1.25) with a diveplane angle (6s) of -20 degrees and rudder angle (6r) of 0 
degrees. The steady state results of the numerical integration method match the linearized 
dynamic results exactly. Vehicle has an out of plane motion in this case. As can be seen 
in Figure 126. the sway velocity V converged to a nonzero steady state value. 



SURGE VELOCITY 




Figure 126. Time History of Surge Velocity for Case 2. 
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SWAY velocity 




Figures 128 and 129 show a plot of surge velocity (U), and sway velocity (V) 
versus time respectively for the Xob=-0.75 % L with a diveplane angle (6s) of -20 
degrees, and rudder angle (6r) of 0 degrees. Here we can see the periodic out of plane 
behavior of the steady state solution. This is what we expected from the eigenvalue 
analysis that we did in the previous section. 

Figures 130 and 131 show a plot of surge velocity (U), and sway velocity (V) 
versus time respectively for x<-b=0 % L with a dive plane angle (6s) of -20 degrees, and 
rudder angle (6r) of 0 degrees. And once again, the steady state results of the numerical 
integration method match the linearized dynamic results exactly. 
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SURGE VELOCITY 




SWAY VELOCITY 




Figure 129. Time History of Sway Velocity for Case 3. 
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SURGE VELOCITY 




SWAY velocity 




Figure 131. Time History of Sway Velocity for Case 4. 
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X. CONCLUSIONS AND RECOMMENDATIONS 



The problem of steady state response and dynamic stability analysis of submarines 
in free positive buoyancy ascent has been studied. The main conclusions of this work can 
be summarized as follows: 

• Steady state motion is not, in general, restricted to the vertical plane, there can 
exist complicated out-of-plane motions. A combination of analytical results for the 
in plane motions and numerical simulations and continuation was used in order to 
identify all possible steady states. 

• For some steady state solutions there can be pitchfork and hysteresis bifurcations. 
In fact, it was shown that the pitchfork is the primary mechanism for generating 
out of plane motions. 

• The response of the vehicle can be very sensitive to system parameters. In such a 
case, any numerical integration results should be viewed with extreme caution. 

• The bifurcation graphs that we generated are very useful in order to understand the 
steady state solution set of the vehicle for any parameter combination. 

• As a recommendation for future research, continuation methods should be used to 
investigate the nature of the Hopf bifurcation, specifically stability of periodic 
solutions. 

• Furthermore, the effect of nonzero propulsion should be investigated. It is possible 
that small nonzero propulsion forces could affect significantly the bifurcations. 
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APPENDIX 



1. BIFPACK 

Biipack is a program package for calculating bifurcations written in 
FORTRAN by R. Seydel [Ref. 15] 

In this work we used Biipack to find the bifurcation characteristics of steady 
state solutions of submersibles. The resulting set of algebraic equations was solved by 
using PACKA which is devoted to system of nonlinear equations 

T (y, X ) =0 

Here y and f(y,X) are vectors with n components. 

To use Biipack we did not bother about how to program our original f(y,X) 
into the enlarged form that is used by Biipack. This larger frame that embraces the 
original n scalar algebraic equations by Biipack is called FRAMEA in Biipack. To solve 
our problem we duplicated this frame and gave it different a name. This became our 
main program and then we entered n scalar formulas fl> ^2 jfsv .., f„ into subroutine 
FCN. A sample of those two programs is included in the following sections. In 
conjunction with subroutine FCN, we also wrote functions which compute the integral 
values in Equations 8, 9, 11, and 12. These are functions F12, FI3, FIS, FI6 in 
subroutine FCN. 
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2. MAIN PROGRAM 



C SUB2.F0R 

PROGRAM MAIN 

IMPUCIT DOUBLE PREOSION (A-H,0-Z) 
DOUBLE PRECISION EPS 
EXTERNAL FCN 



INTERACTIVE WORK WITH THE PACK-A SUBPACKAGE OF 
BIFPACK 

COMMON/DATA/DS,DB,DELB,XGB,ZGB,XB,ZB,DRS,DRB 
OPEN (11, FILE = ’DATAAr,STATUS = ’OLD’) 

OPEN (8 , FILE = ’DIAGRAM’, STATUS = ’OLD’) 

OPEN (3 , FILE = ’START’,STATUS = ’OLD’) 

OPEN (1 , FILE = ’OPTIONS’, STATUS = ’OLD’) 

OPEN (10, FILE =’FORTRAN,MAT’,STATUS = ’OLD’) 

OPEN (20, FILE =’FORTRANl.MAT’,STATUS = ’OLD’) 

OPEN (15, FILE=’DEN.MAT’,STATUS = ’OLD’) 

OPEN (16, FILE= ’DENI .MAT’,STATUS = ’OLD’) 

PROBLEM DEPENDENT VARIABLES: 

EPS = (REL.) ACCURACY, N = NUMBER OF EQUATIONS 



EPS = l.E-5 
N=9 

JACOBI = 1 
DATA 

PI =4.0*ATAN(1.0) 

WHGHT= 12000.0 

L= 17.425 

DB=0.0*PI/180.0 

DELB = 2. 0*WHGHT/ 100.0 

XGB=-1.0-"L/100.0 

ZGB=0.1 

XB=0.0*L/ 100.0 

ZB=0.0 

DRS=0.0-"PI/180.0 



120 



DRB=0.0*PI/180.0 
WRITE (*,*) ’WHATISDRS’ 
READ (*,*) DRS 
DRS=DRS*PI/180.0 
WRITE (*,*) DRS 



1 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 



WRITE (11,*) ’DS AS PARAMETER ZGB = .1,DB IS 2 PERCENT 
OF WEIGHT’ 

WRITE(11,*) ’ = = = = = = = = = = = = = = = = = = = = = = = = =-.= 
INTTLAL SOLUTION USED FOR THIS RUN 



Y(l)=8.6 

Y(2)=0 

Y(3)=-0.25 

Y(4)=0 

Y(5)=0 

Y(6)=0 

Y(7)=0 

Y(8)=0.53 

Y(9)=0 



U 

V 

w 

p 

Q 

R 

pm 

THETA 

PSIDOT 



WRITE (*,*) ’ Previous data in DIAGRAM to be purged?’ 
WRITE (*,*) ’ ENTER 0 FOR PURGING, ENTER 1 FOR 
APPENDING:’ 

READ (*,*) NPURGE 
IF (NPURGE.EQ.O) REWIND 8 
IF (NPURGE.EQ.1) THEN 
DO 23 1=1,9999 

23 READ (8,’()’,END=29) 

29 CONTINUE 

END IF 



CALL INTERA (N,FCN,EPS, JACOBI) 

WRITE (*,*) ’ PRINT-OUTPUT IN FILE "DATAAl"’ 

STOP 

END 



3. SUBROUTINE FCN. 



C EOMl.FOR 

SUBROUTINE FCN (TDUMMY,Y,F) 

IMPUCIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION Y,F,PAR,A,ETACO,ZETACO,PARCO,TDUMMY 
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1 

2 

1 

2 

3 

4 

5 

6 

7 

8 

9 

1 

2 

C 

C 

C 

c 

1 

c 

c 

c 



DOUBLE PRECISION DB, DR, XG, ZG, YG, XPROP, 12, 13, YB, 
KPROP, 15, XB, ZB, BOY, 16 ,NPROP, U, V, W, 

FI2, n3,FI5,FI6 

DOUBLE PREaSION WEIGHT, IX, lY, IZ, IXY, lYZ, IXZ, L, RHO, 

G, AO, CDO, M, XPP, XQQ, XRR, XPR, XUDOT, 

XWQ, XVP, XVR, XQDS, XQDB, XDRDR, XRDR, XVV, 
XWW,XVDR,XWDS,XWDB,XDSDS,XDBDB,XRES, YPDOT, 
YRDOT,YPQ,YQR,YVDOT,YP,YR,YVQ,YWD,YWR,YV, 
YVW, YDR, YDRB, ZVV, ZDS, ZDB, KPDOT, KVW, 
MQDOT, MPP, MPR, ZQDOT, ZPP, ZPR, ZRR, 
ZWDOT,ZQ,ZVP,ZVR,ZW,NRDOT,NPQ,NQR, KRDOT, 
KPQ, KQR, KVDOT, KP, KR, KVQ, KWP, KWR, KV, 

MRR, MWDOT, MQ, MVP, MVR, MW, MVV, MDS, MDB, 
NPDOT,NVDOT, NP, NR, NVQ, NWP, NWR, NV, NVW, 
NDR, NDRB, XRDRB, XRDRS 
PARAMETER (NDIM=12,NDIML=13) 

DIMENSION Y(NDIML),F(NDIML) 

FOLLOWING THREE STATEMENTS ARE FOR CONTINUATION AN 
BRANCH SWITCHING PURPOSES: 

COMMON/DATA/DS,DB,DELB,XGB,ZGB,XB,ZB,DRS,DRB 
COMMON /CONTRO/ ETACO,ZETACO,PARCO,NCO,NlCO,N2CO, 
KCO,IND2CO,NFLAG 
PAR=Y(N1C0) 

F(NlCO) = Y(KCO)-ZETACO*Y(IND2CO)-ETACO 

GEOMETRIC PROPERTIES 

WEIGHT= 12000.0 

IX=1760.0 

lY =9450.0 

IZ= 10700.0 

DCY=0.0 

IYZ=0.0 

IXZ =0.0 

L= 17.425 

RHO = 1.94 

G=32.2 

YG=0.1 

YB=0.0 

A0=2.0 

CDO =0.0057 
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MASS=WHGHT/G 

SURGE HYDRODYNAMIC COEPnCIENTS 

XPP = 7.030E-03*0.5*RHO*L**4 
XQQ =-1.470E-02*0.5’^RHO*L*M 
XRR = 4.010E-03'*‘0.5*RHO*L’^*4 
XPR = 7.640E-04*0.5*RHO*L**4 
XUDOT =-7.580E-03*0.5*RHO*L**3 
XWQ =-1.920E-01*0.5*RHO*L**3 
XVP =-3.240E-03*0.5*RHO*L**3 
XVR = 1.890E-02*0.5*RHO*L**3 
XQDS = 2.610E-02’^0.5'*‘RHO*L**3 
XQDB =-2.600E-03*0.5*RHO*L’^*3 
XRDR =-8.180E-04*0.5*RHO*L*-^3 
XW = 5.290E-02*0.5*RHO*L**2 
XWW = 1.710E-01*0.5*RHO*L**2 
XVDR = 1.730E-03*0.5*RHO*L**2 
XWDS = 4.600E-02*0.5*RHO*L**2 
XWDB = 9.660E-03*0.5*RHO*L**2 
XDSDS =-M60E-02*0.5*RHO*L**2 
XDBDB =-8.070E-03*0.5*RHO*L**2 
XDRDR =-1.010E-02*0.5*RHO*L**2 
XRES = CD0*0.5*RHO*L**2 

XPROP = XRES*ALPHA**2 
XRDRB=XRDR 
XRDRS=XRDR 
XVDRS=XVDR 
XDRB=XVDR 



SWAY HYDRODYNAMIC COEFFICIENTS 

YPDOT = 1.270E-04*0.5*RHO*L*M 
YRDOT = 1.240E-03*0.5*RHO*L**4 
YPQ = 4.125E-03*0.5*RHO*L**4 
YQR =-6.510E-03*0.5*RHO*L**4 
YVDOT =-5.550E-02*0.5*RHO*L’*'*3 
YP = 3.055E-03*0.5*RHO*L-*'*3 
YR = 2.970E-02*0.5*RHO*L**3 
YVQ = 2.360E-02*0.5*RHO*L**3 
YWP = 2.350E-01*0.5*RHO*L**3 
YWR =-1.880E-02*0.5*RHO*L**3 
YV =-9.310E-02*0.5*RHO*L**2 
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YVW = 6.840E-02*0.5’^RHO*L-^’"2 
YDRS =+2.270E-02-"0.5*RHO*L**2 
YDRB =+2.270E-02*0.5*RHO*L’^’^2 
C 

C HEAVE HYDRODYNAMIC COEFFICIENTS 
C 

ZQDOT =-6.810E-03*0.5*RHO’^L*M 
ZPP = 1.270E-04*0.5*RHO’^L**4 
ZPR = 6.670E-03’^0.5*RHO*L*M 
ZRR =-7.350E-03’^0.5’"RHO*L*M 
ZWDOT =-2.430E-01*0.5*RHO*L**3 
ZQ =-1.350E-01*0.5*RHO*L*’"3 
ZVP =-4.810E-02’^0.5’^RHO*L’^’^3 
ZVR = 4.550E-02*0.5*RHO*L**3 
ZW =-3.020E-0H0.5*RHO*L**2 
ZW =-6.840E-02*0.5*RHO’^L’^’^2 
ZDS =-2.270E-02*0.5*RHO*L»*2 
ZDB =-2.270E-02*0.5*RHO-^L*-"2 
C 

C ROLL HYDRODYNAMIC COEFHCIENTS 
C 

KPDOT =-1.010E-03*0.5*RHO*L**5 
KRDOT =-3.370E-05*0.5*RHO*L**5 
KPQ =-6.930E-05*0.5*RHO*L**5 
KQR = 1.680E-02*0.5*RHO*L’"*5 
KVDOT = 1.270E-04*0.5*RHO’^L’^M 
KP =-1.100E-02*0.5*RHO’"L’^*4 
KR =-8.410E-04*0.5*RHO’^L*M 
KVQ =-5.115E-03*0.5*RHO’"L»*4 
KWP =-1.270E-04*0.5>^RHO*L*M 
KWR = 1.390E-02*0.5-^RHO*L-^M 
KV = 3.055E-03*0.5*RHO*L*^*3 
KVW =-1.870E-01*0.5*RHO*L**3 
C 

C PITCH HYDRODYNAMIC COEFnciENTS 
C 

MQDOT =-1.680E-02*0.5’^RHO*L*’"5 
MPP = 5.260E-05-^0.5*RHO*L-^*5 
MPR = 5.040E-03*0.5*RHO*L*»5 
MRR =-2.860E-03*0.5*RHO*L**5 
MWDOT =-6.810E-02*0.5*RHO’"L*M 
MQ =-6.860E-02*0.5*RHO*L*M 
MVP = 1.180E-03*0.5*RHO*L*M 
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MVR = 1.730E-02*0.5*RHO*L**4 
MW = 9.860E-02*0.5*RHO*L**3 
MW =-2.510E-02*0.5*RHO*L**3 
MDS =-1.113E-02’^0.5*RHO’^L’^*3 
MDB = 1.113E-02’^0.5*RHO*L*»3 

YAW HYDRODYNAMIC COEFHCIENTS 

NPDOT =-3.370E-05*0.5*RHO*L’^"‘5 
NRDOT =-3.400E-03’^0.5’^RHO’^L**5 
NPQ =-2.110E-02*0.5*RHO*L*»5 
NQR = 2.750E-03*0.5*RHO*L’^*5 
NVDOT = 1.240E-03*0.5*RHO*L**4 
NP =-8.405E-04*0.5’^RHO*L**4 
NR =-1.640E-02*0.5*RHO-^L**4 
NVQ =-9.990E-03’^0.5*RHO*L’^M 
NWP =-1.750E-02*0.5*RHO*L*M 
NWR = 7.350E-03»0.5*RHO*L**4 
NV =-7.420E-03*0.5*RHO*L**3 
NVW =-2.670E-02*0.5*RHO*L**3 
NDRS =-1.113E-02*0.5’^RHO*L*’^3 
NDRB = + 1.113E-02*0.5*RHO*L**3 

VARIABLE INmUZATION 



RPM=0.0 

BOY = WEIGHT+DELB 
XG=XB+XGB 
ZG=ZB+ZGB 
ALPHA =0.0 



TDUMMY IS DUMMY 

THIS ROUTINE EVALUATES THE PROBLEM DEFINING FUNCTION 
INPUT IS Y(1),...,Y(N), AND PAR=PARAMETER NOW,F(l),.. 

,F(N) OF THE FUNCTION OF RIGHT HAND SIDE OF F(Y,PAR) ARE 
ENTERED: (SIX DEGREE FREEDOM EQUATION OF MOTION OF C 
SUBMARINE) 

U=Y(1) 

V=Y(2) 

W=Y(3) 

P=Y(4) 

Q=Y(5) 
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R=Y(6) 

PHI=Y(7) 

THETA =Y(8) 

PSID0T=Y(9) 

I2=n2(V,W,Q,R) 

I3=n3(V,W,Q,R) 

I5=FI5(V,W,Q,R) 

I6=n6(V,W,Q,R) 

SWAY =12 
HEAVE=I3 
PITCH =15 
YAW =16 

PARAMETER IS DS FOR THIS RUN 



DS=PAR 
SURGE FORCE 

F(l) = MASS*V"R-MASS*W*Q+MASS’^XG'^Q**2+MASS*XG-^R-^*2- 
& MASS*YG*P*Q-MASS*ZG*P*R+XPP-"P’^’"2 +XQQ*Q**2 +XRR*R*-^2 + 
& XPR*P*R+XWQ*W*Q+XVP*V*P+XVR-^V*R+U*Q* 

& (XQDS*DS+XQDB*DB) + 

& U’^R*(XRDRS*DRS+XRDRB*DRB)+XW*V**2+XWW*W^*2+U*V* 
& (XVDRS*DRS+XDRB-^DRB)+U-^W"(XWDS-^DS+XWDB-^DB)+ 

& (XPSDS*DS**2+XDBDB*DB**2+XDRDR*(DRS**2+DRB**2))-"U*’"2- 

& (WEEGHT-BOY) * DSIN(THETA)+XPROP*RPM*RPM-XRES*U*U 

SWAY FORCE 

F(2) = -MASS*U*R-MASS*XG*P*Q+MASS*YG*R**2-MASS*ZG*Q*R+ 
& YPQ*P*Q+YQR*Q*R+YP*U*P+YR*U*R+YVQ*V^Q+YWP*W*P+ 

& YWR*W*R+YV*U*V+YVW*V*W+YDRS*U**2*DRS+YDRB*U**2* 
& DRB+ (WEIGHT-BOY)* DCOS(THETA)*DSIN(PHI)+MASS*W*P+ 

& MASS*YG*P**2+SWAY 

HEAVE FORCE 

F(3) = MASS*U*Q-MASS*V*P-MASS*XG*P*R-MASS*YG*Q*R+ 

& MASS*ZG*P**2+MASS*ZG*Q**2+ZPP*P**2+ZPR*P*R+ 

& ZRR*R**2+ZQ*U*Q+ZVP*V*P+ZVR*V*R+ZW*U*W+ZVV*V**2-I- 
& HEAVE+ U**2*(ZDS*DS+ZDB*DB)+(WEIGHT-B0Y)* 

& DCOS(THETA)* DCOS(PHI) 
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ROLL MOMENT 

F(4) = -IZ*Q*R+IY*Q*R-IXY*P*R+IYZ’^Q**2-IYZ*R*’*'2+IXZ’*'P*Q+ 
& MASS*YG*U*Q-MASS*YG*V*P-MASS*ZG*W*P+KPQ*P’^Q+ 

& KQR’^Q’^R+ KP’^U*P+KR*U*R+KVQ*V’^Q+KWP*W*P+KWR* 

& W*R+KV*U*V+ KVW*V*W+ 

& (YG*WEIGHT-YB*BOY)*DCOS(THETA)*DCOS(PHI)-(ZG*WHGHT- 
& ZB*BOY)*DCOS(THETA)*DSIN(PHI) +MASS*ZG*U*R 

PITCH MOMENT 



F(5) = -IX’^P*R+IZ*P*R+IXY’"Q*R-IYZ’^P*Q-IXZ’"P’"’^2+IXZ’^R*’^2- 
& MASS*XG*U*Q+MASS*XG*V*P+MASS*ZG*V*R-MASS*ZG*W*Q+ 
& MPP*P**2+MPR*P*R+MRR*R**2+MQ*U*Q+MVP*V*'P+MVR*V*R 
& +MW*U*W+MW*V**2+U**2*(MDS*DS+MDB*DB)-(XG*WEIGHT- 
& XB*BOY)*DCOS(THETA)*DCOS(PHI)- 
& (ZG*WEIGHT-ZB*BOY)*DSIN(THETA)+PITCH 

YAW MOMENT 

F(6) = -IY*P*Q+IX*P*Q+IXY*P**2-IXY*Q’^*2+IYZ*P*R-IXZ*Q*R- 
& MASS*XG*U*R+MASS*XG*W*P-MASS*YG*V*R+MASS*YG*W*Q+ 
& NPQ*P*Q+NQR*Q*R+NP’"U*P+NR*U*R+NVQ*V’"Q+NWP*W’"P+ 
& NWR*W*R+NV*U*V+ NVW*V*W+NDRS*U**2*DRS+NDRB* 

& U**2*DRB+ (XG*WEIGHT-XB*BOY)*DCOS(THETA)*DSIN(PHI) + 

& (YG*WEIGHT-YB*BOY)*DSIN(THETA)+YAW 

KINEMATIC RELATIONS 

F(7) = P+PSIDOT*DSIN(THETA) 

F(8) = Q-PSIDOT*DSIN(PHI)*DCOS(THETA) 

F(9) = R-PSIDOT*DCOS(Pffl)*DCOS(THETA) 

RETURN 
END 

4. SAMPLE RUN 

MARS $ run sub2 
WHAT IS DRS 
0 . 

O.OOOOOOOOOOOOOOOOE+OO 

Previous data in DIAGRAM to be purged? 

ENTER 0 FOR PURGING, ENTER 1 FOR APPENDING: 



127 



0 

THIS IS THE MENU YOU CAN CHOOSE FROM: 



0 : ENTERING STARTING DATA 

1 : CONTINUATION: NO RESTRICTIONS 

2 : CONTINUATION: SETTING OF OPTIONS, THEN RUN 

3 : CONTINUATION, FOLLOWING PREVIOUS OPTIONS 
-1 : ENDING THIS SESSION 

ENTER CHOICE of MODUS: 

0 

ENTER COMPONENT Y( 1): 

4.2 

ENTER COMPONENT Y( 2): 

0 

ENTER COMPONENT Y( 3): 

-.7 

ENTER COMPONENT Y( 4): 

0 

ENTER COMPONENT Y( 5): 

0 

ENTER COMPONENT Y( 6): 

0 

ENTER COMPONENT Y( 7): 

0 

ENTER COMPONENT Y( 8): 

.1745 

ENTER COMPONENT Y( 9): 

0 

ENTER PARAMETER: 

.345 

THIS IS THE MENU YOU CAN CHOOSE FROM: 



0 : ENTERING STARTING DATA 

1 : CONTINUATION: NO RESTRICTIONS 

2 : CONTINUATION: SETTING OF OPTIONS, THEN RUN 

3 : CONTINUATION, FOLLOWING PREVIOUS OPTIONS 
-1 : ENDING THIS SESSION 

ENTER CHOICE of MODUS: 

2 

CURRENT STARTING DATA AT PARAMETER: 0.3450000000000000 
THESE ARE THE CURRENT OPTIONS: 

No: 1 , step length : -0.17406E-01 

No: 2 , index of fixed component : -1 
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No: 3 , level of continuation : 0 

No: 4 , maximum number of continuation steps : 200 

No: 5 , step bounds; Predictor;0.100E+00 , Param,:0.100E+00 , rel. change: 0.10 

No: 6 , stability analysis (yes=l, no=0) : 0 

No: 7 , window : -0.35000E+00 .It. parameter .It. 0.35000E+00 

No: 8 , target point: no final value is set. 

ENTER Number OF OPTION YOU WANT TO CHANGE, OR 
ENTER 0 TO RUN THE CONTINUATION, OR 
ENTER -1 TO RETURN TO MENU: 

1 

ENTER INITIAL STEPSIZE (IN REAL): 

-.005 

ENTER Number OF OPTION YOU WANT TO CHANGE, OR 
ENTER 0 TO RUN THE CONTINUATION, OR 
ENTER -1 TO RETURN TO MENU: 

2 

CHOOSE INDEX FOR INITIAL STEP: 

ENTER 0 FOR CONTIN. WITH RESPECT TO PARAMETER, 

ENTER K FOR FIXING K-TH COMPONENT; 

0 

ENTER Number OF OPTION YOU WANT TO CHANGE, OR 
ENTER 0 TO RUN THE CONTINUATION, OR 
ENTER -1 TO RETURN TO MENU: 

3 

CHOOSE AMONG THREE LEVELS OF STEP CONTROL; 

ENTER 2 FOR FREELY VARIABLE STEPS, 

1 FOR VARIABLE STEP BUT FIXED INDEX, 

0 FOR BOTH STEPLENGTH AND INDEX FIXED; 

0 

ENTER Number OF OPTION YOU WANT TO CHANGE, OR 
ENTER 0 TO RUN THE CONTINUATION, OR 
ENTER -1 TO RETURN TO MENU: 

4 

Enter maximum NUMBER OF CONTINUATION STEPS: 

100 

ENTER Number OF OPTION YOU WANT TO CHANGE, OR 
ENTER 0 TO RUN THE CONTINUATION, OR 
ENTER -1 TO RETURN TO MENU: 

-1 

THIS IS THE MENU YOU CAN CHOOSE FROM: 



0 : ENTERING STARTING DATA 

1 : CONTINUATION: NO RESTRICTIONS 
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2 : CONTINUATION: SETTING OF OPTIONS, THEN RUN 

3 : CONTINUATION, FOLLOWING PREVIOUS OPTIONS 
-1 : ENDING THIS SESSION 

ENTER CHOICE of MODUS: 

2 

CURRENT STARTING DATA AT PARAMETER: 0.3450000000000000 
THESE ARE THE CURRENT OPTIONS: 

No: 1 , step length : -0.50000E-02 

No: 2 , index of fixed component : 0 

No: 3 , level of continuation : 0 

No: 4 , maximum number of continuation steps : 100 

No: 5 , step bounds: Predictor: 0.100E+ 00 , Param.:0.100E+00 , rel.change:0.10 

No: 6 , stability analysis (yes=l, no=0) : 0 

No: 7 , window : -0.35000E+00 .It. parameter .It. 0.35000E+00 

No: 8 , target point: no final value is set. 

ENTER Number OF OPTION YOU WANT TO CHANGE, OR 
ENTER 0 TO RUN THE CONTINUATION, OR 
ENTER -1 TO RETURN TO MENU: 

0 
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